!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Common layout of the time integrators
!!
!! \n
!!
!! The generic ode problem is described in the abstract type
!! <tt>c_ode</tt>. <em>This type is abstract so that each application
!! can choose the most convenient layout for the solution state (such
!! as a one-dimensional array, or a multi-dimensional array, or a
!! derived type).</em> The fields of the abstract type are used by the
!! time integrators; notice that a particular time integrator does not
!! necessarily use all the fields.
!!
!! Each time integrator must be initialized once, then it can be
!! called an arbitrary number of times to advance the solution, and
!! finally must be cleaned up. This is necessary especially for
!! multistep methods, which have to keep track of various time levels.
!!
!! Notice that many interfaces specify <tt>intent(inout)</tt> instead
!! of the more natural <tt>intent(out)</tt> because the specifications
!! of \c c_ode will likely include large allocatable components, and
!! it would be inefficient to reallocate them at each time step.
!!
!! \todo Many improvements will become possible while better
!! support for f2003 will become available. In particular, we note the
!! following points:
!! <ul>
!!  <li> work variables of the time integrators should be module
!!  variables of abstract type, which are allocated according
!!  to the chosen problem during the initialization step. This
!!  requires support of the MOLD keyword in the ALLOCATE statement:
!!  see <a
!!  href="http://software.intel.com/en-us/forums/showthread.php?t=84587&o=a&s=lr">this
!!  bug</a>;
!!  <li> the type bound mathematical operators should be introduced
!!  in the abstract class \c c_ode, linking them to the
!!  corresponding deferred procedures;
!! </ul>
!<----------------------------------------------------------------------
module mod_time_integrators_base

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_time_integrators_base_constructor, &
   mod_time_integrators_base_destructor,  &
   mod_time_integrators_base_initialized, &
   c_ode

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members
 
 !> ODE problem
 !!
 !! \note As a general rule, the passed argument is the one which is
 !! set by the procedure.
 !! \note We avoid functions because this would require frequent
 !! allocations/deallocations of the function results, while the
 !! subroutine arguments can be allocated once and then passed around
 !! thanks to the \c inout intent.
 type, abstract :: c_ode
 contains
  !> Evaluate the right hand side \f$f(t,y)\f$
  procedure(i_rhs),    deferred, pass(tnd) :: rhs
  !> Summation of two state variables
  procedure(i_add),    deferred, pass(z)   :: add
  !> Scalar multiplication of a state variable
  procedure(i_stimes), deferred, pass(y)   :: stimes
  !> Combined summation and multiplication \f$z = x + ry\f$
  !!
  !! \note This module provides a default implementation based in the
  !! deferred procedures \c add and \c stimes, which should be
  !! overridden when efficiency is important.
  procedure, pass(z) :: op
  !> Solution of the implicit problem \f$ x - \sigma f(t,x) = b \f$.
  !!
  !! This subroutine is required only for implicit solvers, in which
  !! case it solves the implicit problem for arbitrary
  !! \f$\sigma\in\mathbb{R}\f$ and right hand side \f$b\f$. An
  !! additional argument \c xl is provided to pass a state which can
  !! be used to approximate \f$f\f$ (typically by linearization), thus
  !! solving
  !! \f{displaymath}{
  !!  x - \sigma \tilde{f}_{x_l}(t,x) = b
  !! \f}
  !! where \f$\tilde{f}_{x_l}\f$ approximates \f$f\f$ in a
  !! neighbourhood of \f$x_l\f$.
  procedure, pass(x) :: solve
  !> assignment (should be overloaded for better efficiency)
  procedure, pass(y) :: assign
 end type c_ode

 !> tendency computation
 abstract interface
  subroutine i_rhs(tnd,t,uuu)
   import :: wp, c_ode
   implicit none
   real(wp),     intent(in)    :: t   !< time level
   class(c_ode), intent(in)    :: uuu !< present state
   class(c_ode), intent(inout) :: tnd !< tendency
  end subroutine i_rhs
 end interface

 !> summation
 abstract interface
  pure subroutine i_add(z,x,y)
   import :: c_ode
   implicit none
   class(c_ode), intent(in)    :: x, y
   class(c_ode), intent(inout) :: z !< z = x+y
  end subroutine i_add
 end interface

 !> scalar multiplication
 abstract interface
  pure subroutine i_stimes(y,r,x)
   import :: wp, c_ode
   implicit none
   real(wp),     intent(in)    :: r
   class(c_ode), intent(in)    :: x
   class(c_ode), intent(inout) :: y !< y = r*x
  end subroutine i_stimes
 end interface

 ! private members

! Module variables

 ! public members
 logical, protected ::               &
   mod_time_integrators_base_initialized = .false.
 ! private members
 character(len=*), parameter :: &
   this_mod_name = 'mod_time_integrators_base'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_time_integrators_base_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
       (mod_kinds_initialized.eqv..false.   ) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_time_integrators_base_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_time_integrators_base_initialized = .true.
 end subroutine mod_time_integrators_base_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_time_integrators_base_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_time_integrators_base_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_time_integrators_base_initialized = .false.
 end subroutine mod_time_integrators_base_destructor

!-----------------------------------------------------------------------
 
 subroutine solve(t,x,sigma,b,xl)
  real(wp),     intent(in) :: t, sigma
  class(c_ode), intent(in) :: b, xl
  class(c_ode), intent(inout) :: x
   ! override this subroutine if you use implicit integrators
 end subroutine solve
 
!-----------------------------------------------------------------------

 !> \todo \c tmp should be allocated locally
 pure subroutine op(z,x,r,y,tmp)
  real(wp),     intent(in) :: r
  class(c_ode), intent(in) :: x, y
  class(c_ode), intent(inout) :: z, tmp

!  class(c_ode), allocatable :: tmp

   ! the source keyword allocates any allocatable component
!   allocate(tmp,source=y)
   call tmp%stimes( r , y )
   call z%add( x , tmp )
!   deallocate(tmp)
 end subroutine op
 
!-----------------------------------------------------------------------

 pure subroutine assign(y,x)
  class(c_ode), intent(in)    :: x
  class(c_ode), intent(inout) :: y

   call y%stimes( 1.0_wp , x )

 end subroutine assign
 
!-----------------------------------------------------------------------

end module mod_time_integrators_base

